﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class provides a Cholesky decomposition for symmetric matrices to the form GG^T.
    /// </summary>
    [Serializable]
    public class CholeskyDecomposition : AbstractDecomposition
    {
        /// <summary>
        /// The source matrix to decompress.
        /// </summary>
        private SquareMatrix sourceMatrix;

        /// <summary>
        /// The result of the decomposition.
        /// </summary>
        private SquareMatrix l;

        /// <summary>
        /// Initializes a new instance of the <see cref="CholeskyDecomposition"/> class.
        /// </summary>
        /// <param name="sourceMatrix">The source matrix to decompress.</param>
        public CholeskyDecomposition(SquareMatrix sourceMatrix)
        {
            this.sourceMatrix = sourceMatrix;
        }

        /// <summary>
        /// Gets or sets the source matrix to decompress.
        /// </summary>
        /// <value>The source matrix to decompress.</value>
        public SquareMatrix SourceMatrix
        {
            get { return sourceMatrix; }
            set { sourceMatrix = value; }
        }

        /// <summary>
        /// Gets the result of the decomposition.
        /// </summary>
        /// <value>The result of the decomposition.</value>
        public SquareMatrix L
        {
            get { return l; }
        }

        /// <summary>
        /// Executes the Cholesky decomposition for only symmetric matrices.
        /// </summary>
        public override void ExecuteDecomposition()
        {
            if (!this.sourceMatrix.IsSymmetricMatrix)
            {
                throw new ArgumentException("Matrix has to be symmetric.");
            }

            this.l = new SquareMatrix(this.sourceMatrix.Rows);

            for (int i = 0; i < this.sourceMatrix.Rows; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    double sum = 0.0;

                    for (int k = 0; k < j; k++)
                    {
                        sum += this.L.GetValueAtPosition(i, k) * this.L.GetValueAtPosition(j, k);
                    }

                    if (i == j)
                    {
                        this.L.SetValueAtPosition(i, i,
                                                  Math.Sqrt(this.sourceMatrix.GetValueAtPosition(i, i) - sum));
                    }
                    else
                    {
                        this.L.SetValueAtPosition(i, j, 1.0 / this.l.GetValueAtPosition(j, j)
                                                        * (this.sourceMatrix.GetValueAtPosition(i, j) - sum));
                    }
                }

                if (this.L.GetValueAtPosition(i, i) <= 0)
                {
                    throw new ArithmeticException("Matrix not positive definite.");
                }
            }
        }
    }
}